#include "UniformityAnalysis.h"

UniformityAnalysis::UniformityAnalysis(void)
: _numGroups(0), _width(0), _height(0), _seedlingType(0), _seedlingsFound(0), _maxNumSeedlings(0)
{
    _millimetersPerPixel = 0.0;
    _seedlingBlobRect = cvRect(0, 0, 0, 0);
}

UniformityAnalysis::~UniformityAnalysis(void) {
}

void UniformityAnalysis::testAdaptiveThreshold(int blockSize, int param1, int medianBlockSize, int thresh) {	
    IplImage *binaryImg = cvCreateImage(cvSize(_width, _height), IPL_DEPTH_8U, 1);

    cvNamedWindow("ADAPTIVE THRESH", 0);

    cvAdaptiveThreshold(_rChannel, binaryImg, 255, CV_ADAPTIVE_THRESH_MEAN_C, CV_THRESH_BINARY, blockSize, param1);
    cvSmooth(binaryImg, binaryImg, CV_GAUSSIAN, medianBlockSize, medianBlockSize, 0, 0);
    cvShowImage("ADAPTIVE THRESH", binaryImg);
    cvWaitKey();

    cvThreshold(binaryImg, binaryImg, thresh, 255, CV_THRESH_BINARY);
    cvShowImage("ADAPTIVE THRESH", binaryImg);
    cvWaitKey();

    invertImage(binaryImg, binaryImg);
    cvShowImage("ADAPTIVE THRESH", binaryImg);
    cvWaitKey();

    cvDestroyWindow("ADAPTIVE THRESH");
    cvReleaseImage(&binaryImg);
}

void UniformityAnalysis::initialize(IplImage* img, int numGroups, int type, int maxSeedlings) {
    // Initialize the internal variables
    _width = img->width;
    _height = img->height;
    _seedlingList.clear();
    _numGroups = numGroups;
    _seedlingType = type;
    _maxNumSeedlings = maxSeedlings;
    _currentSegmentId = 1;

    // Create the necessary image channel holders
    _bChannel = cvCreateImage(cvSize(_width, _height), IPL_DEPTH_8U, 1);
    _gChannel = cvCreateImage(cvSize(_width, _height), IPL_DEPTH_8U, 1); 
    _rChannel = cvCreateImage(cvSize(_width, _height), IPL_DEPTH_8U, 1);

    // Split RGB image into separate channels
    cvSplit(img, _bChannel, _gChannel, _rChannel, 0);

    // Initialize the segment id image that will be used in later stages
    _singleChannelWidthStep =_bChannel->widthStep;
    _segmentIdImage = new unsigned int[_height * _singleChannelWidthStep];
}

void UniformityAnalysis::destroy(void) {
    // Release the memory used by each hsv channel image
    cvReleaseImage(&(_bChannel));
    cvReleaseImage(&(_gChannel));
    cvReleaseImage(&(_rChannel));
    cvReleaseImage(&(_seedlings));

    // Clear the seedling groups map
    if(_groups.size() > 0) {
        _groups.clear();
    }

    // Clear the seedling groups key set
    if(_groupKeys.size() > 0) {
        _groupKeys.clear();
    }

    // Release segment id image memory
    delete [] _segmentIdImage;
}

double UniformityAnalysis::processImage(IplImage* img, int numGroups, int type, int maxSeedlings) {
    double uniformity = 0.0;

    // Initialize image channels
    initialize(img, numGroups, type, maxSeedlings);

    // Compute the conversion ratio for millimeters to pixels
    //computeConversionRatio();

    // Measure the seedlings within the image
    extractSeedlings();

    // Compute the uniformity
    //uniformity = computeUniformity();

    // Release internal images
    destroy();

    return uniformity;
}

/* method computeCoversionRatio
* Calculates the pixel to millimeter conversion ratio.
* @return the floating-point ratio of how many pixels per millimeter
*/
void UniformityAnalysis::computeConversionRatio(void) {
    // Reset the ratio
    _millimetersPerPixel = 0.0;

    // Create temporary windows and images for display/storage
    //cvNamedWindow("Color Lines", 0);
    IplImage *bSkel = cvCreateImage(cvSize(_width/4, _height/4), IPL_DEPTH_8U, 1);
    cvSmooth(_bChannel, _bChannel, CV_MEDIAN, 3, 3, 0, 0);
    //cvShowImage("Color Lines", _bChannel);

    // Threshold the blue image
    cvAdaptiveThreshold(_bChannel, _bChannel, 255, CV_ADAPTIVE_THRESH_MEAN_C, CV_THRESH_BINARY_INV, 125, 30);
    cvResize(_bChannel, bSkel, CV_INTER_CUBIC);
    bSkel = skeletonizeZS(bSkel);
    cvShowImage("Color Lines", _bChannel);


    // Detect lines within the blue skeleton image using the Hough transform
    CvMemStorage* storage = cvCreateMemStorage(0);
    int defaultThreshold = 200;
    CvSeq* lines = cvHoughLines2(bSkel, storage, CV_HOUGH_STANDARD, 1, CV_PI/180, defaultThreshold, 0, 0);
    while(lines->total < 2) {
        defaultThreshold -= 10;
        lines = cvHoughLines2(bSkel, storage, CV_HOUGH_STANDARD, 1, CV_PI/180, defaultThreshold, 0, 0);
    }
    IplImage* color_dst = cvCreateImage( cvGetSize(bSkel), 8, 3 );
    std::tr1::unordered_multimap<double, double> lineMap;
    std::vector<double> lineMapKeys;
    int i;
    for(i = 0; i < lines->total; i++) {
        // Get the line
        float* line = (float*)cvGetSeqElem(lines,i);

        // Get its components, rho and theta
        float rho = line[0];
        float theta = line[1];

        // Compute the cosine, a, and sine, b, of theta
        double a = cos(theta), b = sin(theta), m = 0, c = 0;

        // Use a and b to compute the slope and intercept of the line
        m = - a / b;
        c = rho / b;

        // Save the slope and intercept for further processing
        lineMap.insert(std::make_pair(m, c));
        lineMapKeys.push_back(m);

        // Compute coordinates for displaying/drawing the line
        CvPoint pt1, pt2;
        double x0 = a*rho, y0 = b*rho;
        pt1.x = cvRound(x0 + 1000*(-b));
        pt1.y = cvRound(y0 + 1000*(a));
        pt2.x = cvRound(x0 - 1000*(-b));
        pt2.y = cvRound(y0 - 1000*(a));
        cvLine( color_dst, pt1, pt2, CV_RGB(255,0,0), 1, 8 );		
    }
    cvShowImage("Color Lines", color_dst);

    // Find the set containing the most parallel lines
    unsigned int mapKeyIndex = 0, maxSetIndex = 0, curSetCount = 0;
    unsigned int maxSetSize = lineMap.count(lineMapKeys.at(0));
    for(mapKeyIndex = 1; mapKeyIndex < lineMapKeys.size(); mapKeyIndex++) {
        curSetCount = lineMap.count(lineMapKeys.at(mapKeyIndex));
        if(maxSetSize < curSetCount) {
            maxSetIndex = mapKeyIndex;
            maxSetSize = curSetCount;
        }
    }

    // Create a sorted vector containing the candidate lines
    double slope = lineMapKeys.at(maxSetIndex);
    std::pair<std::tr1::unordered_multimap<double, double>::const_iterator,
        std::tr1::unordered_multimap<double, double>::const_iterator> 
        mri = lineMap.equal_range(slope);
    std::vector<double> intercepts;
    for(; mri.first != mri.second; ++mri.first) {
        intercepts.push_back((*mri.first).second);
    }

    // Sort the intercepts as the shortest distances
    // will be the ones closest to each other
    std::sort(intercepts.begin(), intercepts.end());
    unsigned int index;
    std::cout << "Sorted intercepts: " << endl;
    for(index = 0; index < intercepts.size(); index++) {
        std::cout << intercepts.at(index) << endl;
    }

    // Find the closest two lines within this set of parallel lines
    unsigned int interceptIndex;
    double minDistance = 0, curDistance = 0;
    std::pair<unsigned int, unsigned int> minDistIntercepts = std::make_pair(0, 1);
    minDistance = intercepts.at(1) - intercepts.at(0);
    for(interceptIndex = 1; interceptIndex < intercepts.size() - 1; interceptIndex++) {
        curDistance = intercepts.at(interceptIndex + 1) - intercepts.at(interceptIndex);
        if(curDistance < minDistance) {
            minDistIntercepts.first = interceptIndex;
            minDistIntercepts.second = interceptIndex + 1;
            minDistance = curDistance;
        }
    }
    std::cout << "Closest lines are: (" << slope << ", " << intercepts.at(minDistIntercepts.first)
        << ") and (" << slope << ", " << intercepts.at(minDistIntercepts.second) << ")" << endl;

    // Compute the millimeters per pixel ratio
    _millimetersPerPixel = 10.0f / minDistance;
    std::cout << "Millimeters per Pixel: " << _millimetersPerPixel << endl;

    // Release images and windows
    cvReleaseImage(&bSkel);
}


/* method measureSeedlings
* Calculates the length in millimeters of all detected seedlings.
* It then uses these lengths to group the seedlings to the closest 5 mm group.
*/
void UniformityAnalysis::extractSeedlings(void) {
    // Create the image which will hold the seedling skeletons
    IplImage *temp = cvCreateImage(cvSize(_width, _height), IPL_DEPTH_8U, 1);
    cvZero(temp);

    // Threshold the green image and skeletonize it
    cvThreshold(_gChannel, _gChannel, 145, 255, CV_THRESH_BINARY);

    // Save the segments image for viewing
    cvSaveImage("images\\Thresholded-145\\t145.png", _gChannel);

    // Detect seedling blobs within the image to determine seedling location
    //CBlobResult blobs = CBlobResult(gSkel, NULL, 0);
    CBlobResult blobs = CBlobResult(_gChannel, NULL, 0);
    CBlob *currentBlob;

    // Go through the blobs
    for (int i = 0; i < blobs.GetNumBlobs(); i++) {  
        currentBlob = blobs.GetBlob(i);  
        CvRect blobBounds = currentBlob->GetBoundingBox();
        double height = blobBounds.height, width = blobBounds.width;

        // Only consider long, thin blobs as these are most likely seedlings
        if(height > HEIGHT_THRESH && height > MIN_RATIO * width) {

            // Add the blob to the extracted seedlings image for further processing
            currentBlob->FillBlob(temp, cvScalar(255,0,0,0), 0, 0);
        }
    }

    // Save the result	
    _seedlings = skeletonizeZS(temp);
    cvSaveImage("images\\Skeletonized\\skeletons.png", _seedlings);

    // Go through the blobs
    blobs = CBlobResult(_seedlings, NULL, 0);
    for (int i = 0; i < blobs.GetNumBlobs(); i++) {  
        currentBlob = blobs.GetBlob(i);  
        CvRect blobBounds = currentBlob->GetBoundingBox();
        double height = blobBounds.height, width = blobBounds.width;

        // Only consider long, thin blobs as these are most likely seedlings
        if(height > HEIGHT_THRESH && height > MIN_RATIO * width) {
            // Add the blob to a list for further processing
            _seedlingList.push_back(currentBlob);
        }
    }

    // Process the seedling blobs found during extraction
    processBlobs(_seedlingList);

    // Release images and destroy windows
    cvReleaseImage(&temp);
}


void UniformityAnalysis::setPixels(IplImage *img, std::vector<coordinate> &pixels, CvScalar color) {
    // Iterate through all pixels in the vector pixels and set
    // the color for those pixels in image img to color color
    coordinate nextPixel;
    for(unsigned int i = 0; i < pixels.size(); i++) {
        nextPixel = pixels.at(i);
        for(int j = 0; j < img->nChannels; j++) {
            ((unsigned char *)(img->imageData + nextPixel.second * img->widthStep))[img->nChannels * nextPixel.first + j] = (unsigned int) color.val[j];
        }
    }
}


void UniformityAnalysis::addSeedlingsToImage(CompleteSeedlings *completeSeedlings, IplImage *completeSeedlingImage) {
    // Go through the segments and draw them onto the segments image
    for(unsigned int i = 0; i < completeSeedlings->size(); i++) {
        for(unsigned int j = 0; j < completeSeedlings->at(i)->numberOfRoots(); j++) {
            setPixels(completeSeedlingImage, completeSeedlings->at(i)->getRootSegmentAt(j)->getSegment(), colors[i % 20]);
        }
    }
}

void UniformityAnalysis::addToImage(std::vector<SeedlingSegment*> *seedlingSegments, IplImage *segmentsImage) {
    // Go through the segments and draw them onto the segments image
    int numSegments = seedlingSegments->size() - 1;
    stringstream ss;
    for(unsigned int i = 0; i < seedlingSegments->size(); i++) {
        setPixels(segmentsImage, seedlingSegments->at(i)->getSegment(), colors[i % 20]);
    }
}

void UniformityAnalysis::addSegmentsToIdentifierImage(std::vector<SeedlingSegment*> *seedlingSegments) {
    // Pointer to the current segment and its unique id
    unsigned int uniqueId;
    std::vector<coordinate> currentSegment;

    // Go through the segments and draw them onto the segments identifier image
    for(unsigned int i = 0; i < seedlingSegments->size(); i++) {
        // Get the current segment and its unique id
        uniqueId = seedlingSegments->at(i)->_segmentGlobalId;
        currentSegment = seedlingSegments->at(i)->getSegment();

        // Go through the segment component pixel coordinates and add them to the identifier image
        for(unsigned int j = 0; j < currentSegment.size(); j++) {
            _segmentIdImage[currentSegment[j].second * _singleChannelWidthStep + currentSegment[j].first] = uniqueId;
        }
    }
}

void UniformityAnalysis::processBlobs(vector<CBlob*> &seedlings) {
    // Sort the vector of blobs by their starting x position from left to right
    sort(seedlings.begin(), seedlings.end(), sortByX);

    // Pointer to the current blob
    CBlob* currentBlob;

    // Create a temporary color image to display all the detected segments
    IplImage *segmentsImage = cvCreateImage(cvGetSize(_seedlings), IPL_DEPTH_8U, 3);
    IplImage *completeSeedlingsImage = cvCreateImage(cvGetSize(_seedlings), IPL_DEPTH_8U, 3);
    cvSetZero(segmentsImage);
    cvSetZero(completeSeedlingsImage);

    // Compute the width threshold for determining overlap (crude check)
    int singleSeedlingWidth = seedlings[0]->GetBoundingBox().width;
    if(seedlings.size() > 1) {
        singleSeedlingWidth = computeBlobWidthThreshold(seedlings);
    }

    // The number of seedlings in the current blob
    int currentNumSeedlings = 1;
    std::vector<range> segmentStartXRanges;

    // Iterate through seedling skeleton blobs to determine the
    // seedlings that be split up, otherwise process them normally
    for(unsigned int i = 0; i < seedlings.size(); i++) {
        // Get the current blob
        currentBlob = seedlings.at(i);

        // Find all segments and store them into a seedling segment map
        SegmentVectorGraphPair* seedlingSegments = findAllSegments(currentBlob);        

        // Prior to sorting the seedling segment vector by Y coordinates, make sure that all seedlings are growing
        // in the same direction. Thus, if a seedling's starting Y coordinate is less than its ending Y coordinate,
        // then reverse the seedling segment and make the old end coordinate the new start coordinate.
        correctSeedlingEndpoints(seedlingSegments->first);

        // Sort the vector so that the seedlings are sorted in descending order by starting Y coordinate
        std::sort(seedlingSegments->first->begin(), seedlingSegments->first->end(), sortByMaxY);

        // Add them to the segments image for verification purposes
        addToImage(seedlingSegments->first, segmentsImage);

        // Add the segments to the segment identifier image for GUI use
        addSegmentsToIdentifierImage(seedlingSegments->first);

        // Determine if there is overlap using the singleSeedlingWidth as a threshold
        if(currentBlob->GetBoundingBox().width >= singleSeedlingWidth) {
            // Find the number of seedlings by scanning the blob
            segmentStartXRanges = numSeedlingsInBlob();
            currentNumSeedlings = segmentStartXRanges.size();
        }
        else {
            currentNumSeedlings = 1;
        }

        // Console output for verification
        std::cout << "Number of segments in this blob:  " << seedlingSegments->first->size() << std::endl;

        // Print out the expected number of seedlings for this blob
        std::cout << "*******       Expected Seedlings: " << currentNumSeedlings << "     *******" << std::endl;

        // Process the segments to combine them into complete seedlings
        CompleteSeedlings *currentSeedlings = createSeedlings(seedlingSegments, segmentStartXRanges);
        //CompleteSeedlings* currentSeedlings = createSeedlingsFromSegments(seedlingSegments);
        //CompleteSeedlings* currentSeedlings = createSeedlings(seedlingSegments);

        // Console output for verification
        std::cout << "Number of seedlings in this blob: " << currentSeedlings->size() << std::endl << std::endl;
        //printSeedlingInfo(currentSeedlings);

        // Show the complete seedlings image
        addSeedlingsToImage(currentSeedlings, completeSeedlingsImage);

        // Delete the old seedling segments map
        delete seedlingSegments;
    }

    // Save the segments image for viewing
    cvSaveImage("images\\Seedling Segments\\seedling-segments.png", segmentsImage);

    // Save the complete seedlings image for viewing
    cvSaveImage("images\\Complete Seedlings\\complete-seedling.png", completeSeedlingsImage);

    // Release the temporary images and other memory used
    cvReleaseImage(&segmentsImage);
    cvReleaseImage(&completeSeedlingsImage);
}


CompleteSeedlings* UniformityAnalysis::createSeedlings(SegmentVectorGraphPair *svg, std::vector<range> segmentStartXRanges) {
    // Create vector which will contain complete seedlings
    CompleteSeedlings *seedlings = new CompleteSeedlings();

    // Get the expected number of seedlings
    int expectedSeedlings = segmentStartXRanges.size();

    // Get the seedling segment vector and graph from the pair
    SeedlingSegmentVector *segmentVector = svg->first;
    SeedlingSegmentGraph *segmentGraph = svg->second;

    // Create a storage container for the seedling segments that most likely 
    // compose the beginning and ending of complete seedlings in this blob.
    std::vector<std::pair<SeedlingSegment*, SeedlingSegment*>> likelyPairings;

    // Compute the y coordinate threshold for pairing seedling segments
    int yUpperEndThreshold = (_seedlingBlobRect.y + _seedlingBlobRect.height) >> 3,
        yLowerEndThreshold = (_seedlingBlobRect.y + _seedlingBlobRect.height) >> 1;

    // Local pointer variables
    SeedlingSegment *lowerEnd = NULL, *likelyLowerEnd = NULL,
                    *upperEnd = NULL, *likelyUpperEnd = NULL,
                    *nullSeedling = NULL;
    coordinate curCoord;
    range currentRange;

    // Using the start segment x ranges, determine the starting segments
    while(!segmentStartXRanges.empty()) {
        // Grab the current range
        currentRange = segmentStartXRanges.back();

        // Remove it from further consideration
        segmentStartXRanges.pop_back();

        // Go through the segments and find the start segment that best fits into the current range
        for(size_t s = 0; s < segmentVector->size(); s++) {
            // Get the current lower end
            lowerEnd = segmentVector->at(s);

            // Get the min coordinate which will be the start of the segment
            curCoord = lowerEnd->getMinCoordinate();

            // Make sure the segment is from the lower half of the seedling
            // blob and is in between the current acceptable x range 
            if(curCoord.second > yLowerEndThreshold
                && currentRange.first <= curCoord.first
                && curCoord.first <= currentRange.second) {
                    if(likelyLowerEnd == NULL) {
                        likelyLowerEnd = lowerEnd;
                    }
                    else if(likelyLowerEnd->getMinCoordinate().second < curCoord.second) {
                        likelyLowerEnd = lowerEnd;
                    }
            }
        }

        // Add best choice for the current range to the list of likely pairings
        likelyPairings.push_back(std::make_pair(likelyLowerEnd, nullSeedling));

        // Reset the likely lower end pointer
        likelyLowerEnd = NULL;
    }

    // Determine the most likely pairings of start and end segments by considering x distance    
    double distanceMetric = INT_MAX, curDistance = 0;
    for(size_t i = 0; i < likelyPairings.size(); i++) {
        // Get the lower end candidate of the current seedling
        lowerEnd = likelyPairings[i].first;

        // Get the max coordinate endpoint of the lower end
        coordinate lowerEndMaxCoordinate = lowerEnd->getMaxCoordinate();

        // Reset xDistance
        distanceMetric = INT_MAX;

        // reset the likely upper end
        likelyUpperEnd = NULL;

        // Find the most likely upper end of the current seedling's lower end
        for(size_t j = 0; j < segmentVector->size(); j++) {
            // Get the current upper end candidate of the current seedling
            upperEnd = segmentVector->at(j);

            // Skip if upperEnd is the same as lowerEnd
            if(upperEnd == lowerEnd) {
                continue;
            }

            // Get the upper end max coordinate
            coordinate upperEndMaxCoordinate = upperEnd->getMaxCoordinate(),
                       upperEndMinCoordinate = upperEnd->getMinCoordinate();

            // Check to make sure the upper end segment is actually above the lower end segment before doing further processing
            if(lowerEndMaxCoordinate.second <= upperEndMaxCoordinate.second
                && lowerEndMaxCoordinate.second <= upperEndMinCoordinate.second) {
                    continue;
            }

            bool upperMaxJunction = isJunctionPixel(upperEndMaxCoordinate),
                 upperMinJunction = isJunctionPixel(upperEndMinCoordinate);

            // Only attempt to pair it if the segment is in the top part of the seedling blob and is a terminal segment
            // Removed from if statement:  && !isJunctionPixel(upperEnd->getSegment().back()) 
            if(upperEndMaxCoordinate.second < yUpperEndThreshold) {

                // Compute the x distance between the upper and lower ends
                int upperEndMaxCoordinateX = upperEndMaxCoordinate.second == upperEnd->getSegment().back().second ? upperEnd->getSegment().back().first : upperEndMaxCoordinate.first;
                double lowerxDistance = LOWER_X_DISTANCE_RATIO * abs(lowerEndMaxCoordinate.first - upperEndMinCoordinate.first);
                double upperxDistance = UPPER_X_DISTANCE_RATIO * abs(lowerEndMaxCoordinate.first - upperEndMaxCoordinateX);
                double yDistance = Y_DISTANCE_RATIO * abs(lowerEndMaxCoordinate.second - upperEndMaxCoordinate.second);
                curDistance = lowerxDistance + upperxDistance - yDistance;

                // Determine if the current upper end is a better pairing for the current lower end
                if(distanceMetric > curDistance) {
                    distanceMetric = curDistance;
                    likelyUpperEnd = upperEnd;
                }
                else if(distanceMetric == curDistance && likelyUpperEnd->getMaxCoordinate().second > upperEnd->getMaxCoordinate().second) {
                    distanceMetric = curDistance;
                    likelyUpperEnd = upperEnd;
                }
            }
        }

        // If the lower end segment has an inflection point, then it likely is a complete seedling by itself
        if(lowerEnd->hasInflectionPoint() || likelyUpperEnd == NULL) {
            // Update the current likely pairing accordingly
            likelyPairings[i].second = lowerEnd;
        }
        // Otherwise an appropriate endpoint was found for this lower end segment
        else {
            // Update the current likely pairing accordingly
            likelyPairings[i].second = likelyUpperEnd;
        }
    }

    // Use the likely pairings to create complete seedlings
    std::pair<SeedlingSegment*, SeedlingSegment*> currentPair;
    for(size_t i = 0; i < likelyPairings.size(); i++) {
        // Grab the current pair
        currentPair = likelyPairings[i];

        // Build a complete seedling from the lower to the upper end
        seedlings->push_back(buildSeedling(segmentGraph, currentPair.first, currentPair.second));
    }

    // Return the result
    return seedlings;
}


CompleteSeedling* UniformityAnalysis::buildSeedling(SeedlingSegmentGraph *segmentGraph, SeedlingSegment *startSegment, SeedlingSegment *endSegment) {
    // Create the complete seedling to be built
    CompleteSeedling *seedling = new CompleteSeedling();

    // Add the start segment to the seedling
    seedling->addToPrimaryRoot(startSegment);

    /* Use the graph to build a complete seedling from the start segment to the end segment */
    // Create a pointer to use while exploring the graph
    SeedlingSegment *nextSegment = startSegment;
    while(nextSegment != endSegment) {
        // Add the next segment to the primary root axis that grows it in the direction of the end segment
        addNextLikelySegment(segmentGraph, seedling, nextSegment, endSegment);

        // Check to see if any progress was made
        if(nextSegment == seedling->getRootSegmentAt(-1)) {
            std::cout << " *******   ERROR:  CANNOT REACH END SEGMENT DESTINATION!! ******* " << std::endl;
            break;
        }

        // Update the next segment
        nextSegment = seedling->getRootSegmentAt(-1);
    }

    // Return the newly built seedling
    return seedling;
}


void UniformityAnalysis::addNextLikelySegment(SeedlingSegmentGraph *segmentGraph, CompleteSeedling *seedling, SeedlingSegment *currentSegment, SeedlingSegment *endSegment) {
    // Get the upper endpoint of the seedling
    coordinate upperEndpoint = currentSegment->getSegment().back(), edgeEndpoint;

    // Compute the junction key to access the other segments connected to the current segment
    _keyCreator.str(std::string());
    _keyCreator << upperEndpoint.first << "-" << upperEndpoint.second;

    // Get the segments connected to the current segment
    SeedlingSegmentEdges edges = segmentGraph->equal_range(_keyCreator.str());

    // Create a temporary iterator to go through the connected segments
    SeedlingSegmentGraph::const_iterator currentEdge;

    // Create a temp SeedlingSegment pointer to use for comparison
    SeedlingSegment *currentBestSegment = NULL;

    // True if the current edge is terminal, false otherwise
    bool isCurrentEdgeTerminal = false;

    // Iterate through all the segments to find the one that grows the primary root towards the end segment
    for(currentEdge = edges.first; currentEdge != edges.second; ++currentEdge) {
        // Skip the segment connecting to this junction that is already part of the seedling
        if(seedling->isPartOfSeedling(currentEdge->second)) {
            continue;
        }

        // Automatically terminate when the current edge is the end segment
        if(currentEdge->second == endSegment) {
            // Set current best segment to the end segment and exit the loop
            currentBestSegment = endSegment;
            break;
        }

        // Grab the end point for the current edge that is not connected to the seedling primary axis
        edgeEndpoint = currentEdge->second->getSegment().back();
        if(edgeEndpoint.first == upperEndpoint.first && edgeEndpoint.second == upperEndpoint.second) {
            edgeEndpoint = currentEdge->second->getSegment()[0];
        }

        // Determine if the edge is terminal
        isCurrentEdgeTerminal = !isConnectingSegment(segmentGraph, currentEdge->second, _keyCreator.str());

        // If the segment grows upward, check whether it grows the primary root axis towards the end segment
        if(edgeEndpoint.second <= upperEndpoint.second && !isCurrentEdgeTerminal) {
            // Take this segment to be the current best if there isn't already a current best
            if(currentBestSegment == NULL) {
                currentBestSegment = currentEdge->second;
            }
            // Otherwise determine which segment is the better choice by looking at the distance to the end segment upper endpoint
            else {
                // Get the upper coordinates of the current best, edge, and end segments
                coordinate bestSegmentUpperEndpoint = currentBestSegment->getSegment().back(),
                           edgeSegmentUpperEndpoint = currentEdge->second->getSegment().back(),
                           endSegmentUpperEndpoint = endSegment->getSegment().back();

                // Compute the x and y coordinate differences between: a) best and end, and b) edge and end
                int xDiffBest = bestSegmentUpperEndpoint.first - endSegmentUpperEndpoint.first, yDiffBest = bestSegmentUpperEndpoint.second - endSegmentUpperEndpoint.second,
                    xDiffEdge = edgeSegmentUpperEndpoint.first - endSegmentUpperEndpoint.first, yDiffEdge = edgeSegmentUpperEndpoint.second - endSegmentUpperEndpoint.second,
                    bestDistance = xDiffBest * xDiffBest + yDiffBest * yDiffBest,
                    edgeDistance = xDiffEdge * xDiffEdge + yDiffEdge * yDiffEdge;

                // Determine whether the best segment needs to be updated
                if(bestDistance > edgeDistance) {
                    currentBestSegment = currentEdge->second;
                }
            }
        }
        // If the segment grows downward, add it as a secondary root as this cannot be part of the main axis
        // However, there should be restrictions on the size of secondary segments
        if(currentEdge->second->getSegment().size() < (seedling->getRootLength() >> 1)) {
            if(edgeEndpoint.second > upperEndpoint.second || isCurrentEdgeTerminal) {
                seedling->addAsSecondaryRoot(currentEdge->second);
            }
        }
    }

    // Add the best segment to the primary axis of the seedling
    seedling->addToPrimaryRoot(currentBestSegment);

    // Clear the key creator stream
    _keyCreator.str(std::string());
}


/* method numSeedlingsInBlob
* Compute the number of seedlings found in the current blob 
* 
* @return the number of seedlings in the blob currently being processed 
*/
std::vector<range> UniformityAnalysis::numSeedlingsInBlob() {
    // Local variables needed to keep track of row specific
    // seedling counts as well as the scan step between rows
    size_t numSeedlingsInBlob = 1, maxCount = 0, minSeparation, 
        seedlingsInCurrentRow,numBins = _maxNumSeedlings + 1;
    unsigned int curSegmentId = 0, yHalfPoint = (_seedlingBlobRect.y + _seedlingBlobRect.height) >> 1;
    bool separated = true;

    // Use a vector keep track of the x coordinate crossings where white pixels are found
    std::vector<int> xCrossings;

    // Use bins to determine the number seedlings by counting white pixels for each scan line
    unsigned int *whitePixelHist = new unsigned int[numBins];

    // Zero out the histogram array
    memset(whitePixelHist, 0, sizeof(int) * numBins);

    // Declare the pointer to the current row in the image
    unsigned char *currentRow;

    // Begin iterating through the image where this blob is found
    for(unsigned int y = _seedlingBlobRect.y; y < (unsigned) (_seedlingBlobRect.y + _seedlingBlobRect.height); y += Y_SCAN_STEP) {
        // Keep a local x coordinate list for the current row
        std::vector<int> currentRowCrossings;

        // Compute the pointer to the current row
        currentRow = (unsigned char*)_seedlings->imageData + y * _seedlings->widthStep;

        // Reset row specific variables
        separated = true;
        curSegmentId = 0;
        minSeparation = 5;
        seedlingsInCurrentRow = 0;

        // Iterate through the current row 
        for(unsigned int x = _seedlingBlobRect.x; x < (unsigned) (_seedlingBlobRect.x + _seedlingBlobRect.width); x++) {
            // Check whether the current pixel is white
            if(currentRow[x] == 255) {
                // For segments that curve or loop we only want to count the crossing once
                if(curSegmentId == 0) {
                    curSegmentId = _segmentIdImage[y * _singleChannelWidthStep + x];
                }
                // This segment has been encountered before
                else if(curSegmentId == _segmentIdImage[y * _singleChannelWidthStep + x]) {
                    separated = false;
                }
                // As long as this segment has not been encountered before and is separated
                // by at least 5 pixels from the previously encountered white pixel, add it
                if(separated) {
                    seedlingsInCurrentRow++;

                    // Add the x coordinates if found in the second half of the blob rectangle area
                    if(y >= yHalfPoint) {
                        xCrossings.push_back(x);
                    }
                            
                }
            }
            // If it isn't white and a white pixel has already been encountered
            else if(currentRow[x] != 255 && !separated) {
                minSeparation--;
                if(minSeparation == 0) {
                    separated = true;
                    minSeparation = 10;
                }
            }
        }

        // Increment the appropriate bin
        seedlingsInCurrentRow = seedlingsInCurrentRow > (unsigned) _maxNumSeedlings ? _maxNumSeedlings : seedlingsInCurrentRow;
        whitePixelHist[seedlingsInCurrentRow]++;
    }

    // Determine the most likely number of seedlings by finding the maximum in the array
    for(size_t n = 1; n < (unsigned) _maxNumSeedlings + 1; n++) {
        if(maxCount < whitePixelHist[n]) {
            maxCount = whitePixelHist[n];
            numSeedlingsInBlob = n;
        }
    }

    // Declare vector container for x crossing ranges
    std::vector<range> xRanges;

    /* Determine the range of x crossings to better pick the lower end of each seedling segment */
    // First sort the x crossings vector
    std::sort(xCrossings.begin(), xCrossings.end());

    // Iterate through the newly sorted array to figure out the beginning and ending of each block
    size_t blockStart = 0;
    for(size_t i = 0; i < xCrossings.size() - 2; i++) {

        // End of current block found
        if((xCrossings[i + 1] - xCrossings[i]) > MIN_BLOCK_SEP) {                      
            // Add information to xRanges vector
            xRanges.push_back(std::make_pair(
                xCrossings[blockStart] - X_RANGE_RELAXATION,
                xCrossings[i] + X_RANGE_RELAXATION));

            // Update the starting position
            blockStart = i + 1;
        }
    }

    // The last block will most likely not be added in the previous for loop, so add it now
    xRanges.push_back(
        std::make_pair(
            xCrossings[blockStart] - X_RANGE_RELAXATION,
            xCrossings[xCrossings.size() - 1] + X_RANGE_RELAXATION));

    // Check to see if any of the ranges can be merged
    std::vector<range> mergedRanges;
    while((xRanges.size() + mergedRanges.size()) > numSeedlingsInBlob) {
        // Grab the last range in the vector
        range last = xRanges.back();
        xRanges.pop_back();

        // Grab the second to last range in the vector
        range secondToLast = xRanges.back();
        xRanges.pop_back();

        // Check to see if the end point of second to last range falls in the last range
        if(last.first <= secondToLast.second && secondToLast.second <= last.second) {
            xRanges.push_back(std::make_pair(secondToLast.first, last.second));
        }
        else {
            mergedRanges.push_back(last);
            xRanges.push_back(secondToLast);
        }

        // If there is only one range left in xRanges, then
        // no more ranges can be merged. So exit the loop.
        if(xRanges.size() == 1) {
            break;
        }
    }

    // Now add back the merged ranges back to the xRanges vector
    xRanges.insert(xRanges.end(), mergedRanges.begin(), mergedRanges.end());

    // Free local memory used
    delete [] whitePixelHist;

    // Return the result
    return xRanges;
}


/* Compute the width threshold used to determine whether multiple seedlings may be in the current blob 
* 
* @param  blobs  the vector of blobs for the current image being processed
* @return integer value equal to the mean + standard devation of the input blobs widths
*/
int UniformityAnalysis::computeBlobWidthThreshold(std::vector<CBlob*> &blobs) {
    // Local storage for computing the mean and variance without first computing the mean
    double *mean = new double[blobs.size()], *variance = new double[blobs.size()];

    // Local storage for the width vector, used for the median
    std::vector<int> blobWidths;

    // Use memset to zero out the arrays
    memset(mean, 0, sizeof(double) * blobs.size());
    memset(variance, 0, sizeof(double) * blobs.size());

    // Initialize the first value of the variance array prior to starting
    mean[0] = blobs.at(0)->GetBoundingBox().width;

    // Now iterate through the remaining blobs computing the mean and variance incrementally. This approach can be
    // found in Donald Knuth's "The Art of Computer Programming, Vol 2: Seminumerical Algorithms" section 4.2.2.
    int curBlobWidth = 0;
    for(size_t i = 1; i < blobs.size(); i++) {
        curBlobWidth = blobs.at(i)->GetBoundingBox().width;
        blobWidths.push_back(curBlobWidth);
        mean[i] = mean[i - 1] + ((curBlobWidth - mean[i - 1]) / (i + 1));
        variance[i] = variance[i - 1] + ((curBlobWidth - mean[i-1]) * (curBlobWidth - mean[i]));
    }

    // Compute the width threshold using the median and standard deviation
    std::sort(blobWidths.begin(), blobWidths.end());
    double stdDev = std::sqrt(variance[blobs.size() - 1] / (blobs.size() - 1)) / 2;
    double widthThreshold = (blobWidths[blobWidths.size() >> 1] + stdDev + 0.5);

    // Free local memory
    delete [] mean;
    delete [] variance;

    // Return the result
    return (int) widthThreshold;
}



CompleteSeedlings* UniformityAnalysis::createSeedlingsFromSegments(SegmentVectorGraphPair* separatedSegments) {
    // Initialize vector of complete seedlings
    CompleteSeedlings *completeSeedlings = new CompleteSeedlings();

    // Split up the input pair for easier use
    SeedlingSegmentVector *segmentVector = separatedSegments->first;
    SeedlingSegmentGraph *segmentGraph = separatedSegments->second;

    // Prior to sorting the seedling segment vector by Y coordinates, make sure that all seedlings are growing
    // in the same direction. Thus, if a seedling's starting Y coordinate is less than its ending Y coordinate,
    // then reverse the seedling segment and make the old end coordinate the new start coordinate.
    correctSeedlingEndpoints(segmentVector);

    // Prune short terminal edges from the graph
    //pruneShortTerminalEdges(segmentGraph);

    // Sort the vector so that the seedlings are sorted in descending order by starting Y coordinate
    std::sort(segmentVector->begin(), segmentVector->end(), sortByMaxY);

    // Go through all the segments and build complete seedlings
    while(!segmentVector->empty()) {
        if(!isPartOfExistingSeedlings(segmentVector->back(), completeSeedlings)) {
            completeSeedlings->push_back(buildSeedling(segmentGraph, segmentVector->back()));
        }
        segmentVector->pop_back();
    }

    // Return all complete seedlings
    return completeSeedlings;
}

void UniformityAnalysis::correctSeedlingEndpoints(SeedlingSegmentVector* segmentVector) {
    SeedlingSegment* currentSegment;
    for(unsigned int i = 0; i < segmentVector->size(); i++) {
        currentSegment = segmentVector->at(i);

        // Remove duplicates
        for(size_t j = i + 1; j < segmentVector->size(); j++) {
            if(currentSegment == segmentVector->at(j)) {
                segmentVector->erase(segmentVector->begin() + j);
            }
            else {
                break;
            }
        }

        if(currentSegment->getSegment().front().second < currentSegment->getSegment().back().second) {
            std::reverse(currentSegment->getSegment().begin(), currentSegment->getSegment().end());
            currentSegment->setStartCoordinate(currentSegment->getSegment().front());
        }

        currentSegment->checkForInflection();
        currentSegment->updateEndDirection();
        currentSegment->updateStartDirection();
        currentSegment->updateOverallDirection();
    }
}

bool UniformityAnalysis::isPartOfExistingSeedlings(SeedlingSegment *segment, CompleteSeedlings *seedlings) {
    // Assume the segment hasn't been used yet
    bool isPartOfSeedling = false;

    // Look through all seedlings' secondary roots  to 
    // ensure the current segment has not been used yet
    for(unsigned int i = 0; i < seedlings->size(); i++) {
        if(seedlings->at(i)->isPartOfSeedling(segment)) {
            isPartOfSeedling = true;
            break;
        }
    }

    // Return the result
    return isPartOfSeedling;
}

CompleteSeedling* UniformityAnalysis::buildSeedling(SeedlingSegmentGraph* segmentGraph, SeedlingSegment* startSegment) {
    // Create the completed seedling object and add the start segment to it
    CompleteSeedling* seedling = new CompleteSeedling();
    seedling->addToPrimaryRoot(startSegment);

    // Pointer to the next segment to examine
    SeedlingSegment* mainSegment = startSegment;

    // Keep a list of the junctions that should be/have been explored while building the seedling
    std::list<coordinate> junctions;
    std::vector<coordinate> visitedJunctions;

    // Add the first junction to be examined
    junctions.push_back(mainSegment->getSegment().back());

    // Build a vector of all segments connected to this junction, if they exist
    while(!junctions.empty()) {
        // Update the main segment
        mainSegment = seedling->getRootSegmentAt(-1);

        // Get the next junction to check
        coordinate nextJunction = junctions.back();
        junctions.pop_back();
        visitedJunctions.push_back(nextJunction);

        // Build the key string
        _keyCreator << nextJunction.first << "-" << nextJunction.second;

        // Check all connected segments and grow the seedling with the best one
        checkConnectedSegments(segmentGraph, seedling, _keyCreator.str());

        // Clear the key creator string stream
        _keyCreator.str("");


        // Update the main segment and add the next junction pixel that should be checked, if any
        mainSegment = seedling->getRootSegmentAt(-1);
        coordinate endOfMainSegment = mainSegment->getSegment().back();
        bool addJunction = true;

        // Search through the visited junctions to ensure we don't look at the same junction twice
        for(unsigned int i = 0; i < visitedJunctions.size(); i++) {
            if( visitedJunctions[i].first == endOfMainSegment.first &&
                visitedJunctions[i].second == endOfMainSegment.second) {
                    addJunction = false;
                    break;
            }
        }

        // Add the junction if we havent checked it already
        if(addJunction) {
            junctions.push_back(mainSegment->getSegment().back());
        }
    }

    return seedling;
}

void UniformityAnalysis::checkConnectedSegments(SeedlingSegmentGraph* segmentGraph, CompleteSeedling* seedling, std::string &mainSegmentKey) {		
    // Graph iterator used to go through all seedling segment edges
    SeedlingSegmentGraph::const_iterator currentSegment;

    // Get the main segment from the connectedSegments vector
    SeedlingSegment* mainSegment = seedling->getRootSegmentAt(-1);

    // A pair of iterators used to obtain pointers to the seedling
    // segments connected to the current junction key
    SeedlingSegmentEdges edges = segmentGraph->equal_range(mainSegmentKey);

    // Seedling Segment Vector which contains connected segments that are potential candidates
    SeedlingSegmentVector connectedSegments;

    // Go through all the segments connected to the current junction
    for(currentSegment = edges.first; currentSegment != edges.second; ++currentSegment) {
        // Skip the segment connectiong to this junction that is already part of the seedling
        if(seedling->isPartOfSeedling(currentSegment->second)) {
            continue;
        }
        // Consider the segment if it increases the length of the main segment
        if(segmentIncreasesLength(mainSegment, currentSegment->second)) {
            connectedSegments.push_back(currentSegment->second);
        }
        // This segment is most likely a secondary root segment. However, before
        // adding it as a secondary segment, check that it is a terminal segment.
        else if(!isConnectingSegment(segmentGraph, currentSegment->second, mainSegmentKey)){
            seedling->addAsSecondaryRoot(currentSegment->second);
        }
    }

    // If no segments passed the first set of tests, there may be a connecting segment that should be used
    if(connectedSegments.empty()) {
        for(currentSegment = edges.first; currentSegment != edges.second; ++currentSegment) {
            // Consider the segment if it connects the main segment to other segments
            if(isConnectingSegment(segmentGraph, currentSegment->second, mainSegmentKey)) {
                connectedSegments.push_back(currentSegment->second);
            }
        }
    }

    // Add the most likely segment to the main segment based on the angle of the connection
    if(connectedSegments.size() > 1) {
        seedling->addToPrimaryRoot(selectBestSegment(seedling, connectedSegments));
    }
    // If there is only one segment just add it
    else if(!connectedSegments.empty()) {
        seedling->addToPrimaryRoot(connectedSegments[0]);
    }
}

bool UniformityAnalysis::isConnectingSegment(SeedlingSegmentGraph* segmentGraph, SeedlingSegment* testSegment, std::string &mainSegmentKey) {
    // Flag that indicates whether the test segment
    // connects the main segment to other segments
    bool connectingSegment = false;

    // Clear the string stream object used to build the junction keys needed for the graph
    _keyCreator.str(std::string());

    // Build the key for the end point that isn't connected to the main segment
    coordinate* endPoint = &(testSegment->getSegment().front());
    _keyCreator << endPoint->first << "-" << endPoint->second;

    // Check that it is not the main segment key
    if(mainSegmentKey == _keyCreator.str()) {
        // Reset and rebuild the key using the other end point of the test segment
        _keyCreator.str(std::string());
        endPoint = &(testSegment->getSegment().back());
        _keyCreator << endPoint->first << "-" << endPoint->second;
    }

    // Check whether the end points of testSegment both connect to other segments
    SeedlingSegmentEdges edges = segmentGraph->equal_range(_keyCreator.str());

    // Graph iterator used to go through all seedling segment edges
    SeedlingSegmentGraph::const_iterator currentSegment;

    // Go through all the segments connected to the current junction
    int numberOfConnectedSegments = 0;
    for(currentSegment = edges.first; currentSegment != edges.second; ++currentSegment) {
        if(currentSegment->second != testSegment) {
            numberOfConnectedSegments++;
        }
    }

    // If the number of connected segments is > 0, then this is a connecting segment
    if(numberOfConnectedSegments > 0) {
        connectingSegment = true;
    }

    // Return the result
    return connectingSegment;
}


/* method findAllSegments
* Given a seedling blob containing multiple overlapping seedlings,
* this method splits the blob into segments in order to create a
* graph where the joining or junction pixels are vertices and the
* seedling segments are edges.  This graph, known as the Seedling
* Segment Graph is a multimap where the junction pixel coordinate
* in its string form is the key and seedling segments are values.
*
* @param  blob the input CBlob seedling blob
* @return a seedling segment graph used to represent the input blob
*/
SegmentVectorGraphPair* UniformityAnalysis::findAllSegments(CBlob* blob) {
    // Initialize the Seedling Segment Graph, Vector and  Vect/Graph Pair
    // The SSG is a multimap that uses junction pixel coordinates as keys
    SeedlingSegmentGraph *segmentGraph = new SeedlingSegmentGraph(); 
    SeedlingSegmentVector *allSegments = new SeedlingSegmentVector();
    SegmentVectorGraphPair *svgp = new SegmentVectorGraphPair(allSegments, segmentGraph);

    // Initialize an array used to keep track of 
    // visited pixels in the connected component
    initVisitedPixelsArray(blob);


    // Find the start coordinate of the first foreground pixel
    coordinate *startCoordinate = new coordinate(findSeedlingStartPoint(blob)),
        *nextJunction = new coordinate(-1, -1);
    stringstream keyStream;

    // Local space for the junction pixel queue used to build the connected component
    std::list<coordinate*> junctionQueue;

    // Add the first pixel
    junctionQueue.push_back(startCoordinate);

    while(!junctionQueue.empty()) {
        // Create a seedling segment used to accumulate the current segment
        SeedlingSegment* seedling = new SeedlingSegment();

        // Get the next pixel in the junction queue
        coordinate *pixelCoord = junctionQueue.front();

        // Remove it from the queue
        junctionQueue.pop_front();

        // Mark the pixel as visited
        setVisited(pixelCoord->first, pixelCoord->second);

        // Add the new segment to the segment graph
        keyStream << pixelCoord->first << "-" << pixelCoord->second;
        segmentGraph->insert(std::make_pair(keyStream.str(), seedling));
        allSegments->push_back(seedling);
        keyStream.str(std::string());

        // Build the new segment until either a terminal or junction pixel is encountered
        buildSegment(seedling, pixelCoord->first, pixelCoord->second, nextJunction);

        // Assign it an integer unique identifier
        seedling->_segmentGlobalId = _currentSegmentId;

        // Add it to the segment identifier map for quick GUI access support
        _segmentIdentifierMap.insert(std::make_pair(_currentSegmentId, seedling));

        // Increment the segment unique identifier
        _currentSegmentId++;

        // Determine whether a junction was encountered
        if(nextJunction->first != -1 && nextJunction->second != -1) {
            // If one was encountered, add it to the queue for processing.
            // Since a junction has multiple branches, add it to the queue
            // as many times as it has explorable neighbors.
            std::vector<coordinate> neighbors;
            getNeighbors(nextJunction->first, nextJunction->second, neighbors);
            for(unsigned int i = 0; i < neighbors.size(); i++) {
                junctionQueue.push_back(new coordinate(nextJunction->first, nextJunction->second));
            }

            // Add the current seedling segment to the graph under
            // this junction pixel if it hasn't been added already
            keyStream << nextJunction->first << "-" << nextJunction->second;
            if(!graphContainsSegment(segmentGraph, keyStream.str(), seedling)) {
                segmentGraph->insert(std::make_pair(keyStream.str(), seedling));
                allSegments->push_back(seedling);
            }
            keyStream.str(std::string());

            // Reset nextJunction
            nextJunction->first = -1;
            nextJunction->second = -1;
        }
    }

    // Delete the visited pixels scratch space and other dynamically allocated space
    delete [] _visitedPixels;
    delete nextJunction;

    // Return the seedling segment graph & vector pair
    return svgp;
}

/* method buildSegment
* Builds the connected component of the primary axis for the input seedling
* segment looking at the center position (x, y) within the _seedling image.
* 
* @param segment  the seedling segment used to accumulate the connected component
* @param x  the x coordinate of the starting center pixel
* @param y  the y coordinate of the starting center pixel 
* @param junction  the coordinate of the first junction pixel 
*                  encountered while building the segment
*/
void UniformityAnalysis::buildSegment(SeedlingSegment* segment, int x, int y, coordinate *junction) {	

    // Local space for the pixel neighborhood and queue used to build the connected component
    std::list<coordinate> pixelQueue;
    std::vector<coordinate> neighbors;

    // Initialize the local queue
    pixelQueue.push_back(std::make_pair(x, y));

    // While the component queue still has pixels to process
    while(!pixelQueue.empty()) {
        // Get the next pixel to process
        coordinate nextPixel = pixelQueue.front();
        pixelQueue.pop_front();

        // Clear any previous neighbors if any
        neighbors.clear();

        // Get the unvisited connected neighbors of this pixel
        getNeighbors(nextPixel.first, nextPixel.second, neighbors);

        // Grow this segment
        segment->addToSegment(nextPixel.first, nextPixel.second);

        // nextPixel is a junction pixel (e.g. a branching point)
        if(neighbors.size() > 1 && isNotVisited(nextPixel.first, nextPixel.second)) {
            // Check to see if the current pixel is a junction pixel
            if(isJunctionPixel(nextPixel)) {
                // Set the junction's coordinates to this pixel's coordinates
                junction->first = nextPixel.first;
                junction->second = nextPixel.second;
            }
            // Otherwise determine which neighbor is the junction pixel and explore that next
            else {
                // Search for the junction pixel in nextPixel's neighborhood
                for(unsigned int i = 0; i < neighbors.size(); i++) {
                    coordinate neighbor = neighbors.at(i);
                    if(isJunctionPixel(neighbor) && isNotVisited(neighbor.first, neighbor.second)) {
                        pixelQueue.push_back(neighbor);
                        break;
                    }
                }

                // If the search did not yeild a junction pixel, then the neighboring 
                // junction pixel has already been explored. The neighbor that should
                // be expanded next therefore needs to be the one that is connected to
                // the current pixel, but not the junction pixel.
                if(pixelQueue.empty()) {
                    for(unsigned int i = 0; i < neighbors.size(); i++) {
                        coordinate neighbor = neighbors.at(i);
                        if(!isConnectedToJunction(neighbor)) {
                            pixelQueue.push_back(neighbor);
                        }
                    }
                }
            }
        }
        // The pixel being explored is a junction pixel that has been
        // visited already or a pixel with only 1 unvisited neighbor
        else if((!isNotVisited(nextPixel.first, nextPixel.second) && neighbors.size() > 1) || neighbors.size() == 1) {
            // Expand and start exploring its first neighbor
            pixelQueue.push_back(neighbors.at(0));
        }

        // Mark it as visited
        setVisited(nextPixel.first, nextPixel.second);
    }
}

/* method isJunctionPixel
* Computes the junction code for the neighborhood of the input center pixel coordinate.
* If the compute code is one of the junction codes present in the junction sum set, then
* the input center pixel is a junction pixel and isJunctionPixel returns true.
* 
* @param center an std::pair object with the first element being the x coordinate and
*               the second element being the y coordinate.
* @return true if the pixel in question is a junction pixel, false otherwise
*/
bool UniformityAnalysis::isJunctionPixel(coordinate center) {
    // Assume pixel in question is not a junction, as the majority are not
    bool isJunction = false;

    // Initialize the neighborhood junction code sum to zero
    int junctionCode = 0;

    // Compute some coordinates for convenience
    int xMax = _seedlingBlobRect.x + _seedlingBlobRect.width,
        yMax = _seedlingBlobRect.y + _seedlingBlobRect.height,
        xPrev = center.first - 1,
        yPrev = center.second - 1,
        xNext = center.first + 1,
        yNext = center.second + 1;

    // Compute pointers to rows in the _seedlings image for convenience
    unsigned char *currentRow = (unsigned char*)_seedlings->imageData + center.second * _seedlings->widthStep,
        *previousRow = (unsigned char*)_seedlings->imageData + yPrev * _seedlings->widthStep,
        *nextRow = (unsigned char*)_seedlings->imageData + yNext * _seedlings->widthStep;

    /* Compute the sum of pixels in the neighborhood of center */

    // Check Bottom Row of Pixels
    if(yNext < yMax) {
        if(xPrev > 0 && nextRow[xPrev] > 0) {
            junctionCode += 1;
        }
        if(nextRow[center.first] > 0) {
            junctionCode += 2;
        }
        if(xNext < xMax && nextRow[xNext] > 0) {
            junctionCode += 4;
        }
    }

    // Check Current Row of Pixels
    if(xPrev > 0 && currentRow[xPrev] > 0) {
        junctionCode += 128;
    }
    if(xNext < xMax && currentRow[xNext] > 0) {
        junctionCode += 8;
    }

    // Check Top Row of Pixels
    if(yPrev > 0) {
        if(xPrev > 0 && previousRow[xPrev] > 0) {
            junctionCode += 64;
        }
        if(previousRow[center.first] > 0) {
            junctionCode += 32;
        }
        if(xNext < xMax && previousRow[xNext] > 0) {
            junctionCode += 16;
        }
    }

    // Check to see if the computed junction code is a member of the valid junction sums set.
    // If so, then the center pixel is a junction. Otherwise, the center pixel is not.
    if(junctionSumSet.find(junctionCode) != junctionSumSet.end()) {
        isJunction = true;
    }

    return isJunction;
}

/* method isConnectToJunction
* Searches the neighborhood of pixel coordinate center for a junction pixel.
* Returns true if one is found or false if one is not.
* 
* @param center  coordinate of the pixel in question
* @return true   if connected to junction pixel
*         false  otherwise
*/
bool UniformityAnalysis::isConnectedToJunction(coordinate center) {
    bool isConnected = false;

    // Search the neighborhood of pixel coordinate center for a junction pixel
    for(int i = -1; i < 2 && !isConnected; i++) {
        for(int j = -1; j < 2 && !isConnected; j++) {
            isConnected = isJunctionPixel(coordinate(center.first + j, center.second + i));
        }
    }

    return isConnected;
}

/* method graphContainsSegment
* Given a SeedlingSegmentGraph, junction coordinate key, and SeedlingSegment,
* this function checks to see whether the input seedling already exists within
* the input graph.
* 
* @param segmentGraph  the multimap representing the graph of the current seedling blob being processed
* @param junctionKey   the junction coordinate key string used to get the list of seedling segments of
*                      the seedling segment graph connected to the particular junction in question
* @param seedling      the seedling segment we are searching for
* @return  true   if the seedling is in the graph at the location in question, 
*          false  otherwise
*/
bool UniformityAnalysis::graphContainsSegment(SeedlingSegmentGraph *segmentGraph, std::string junctionKey, SeedlingSegment *seedling) {
    bool segmentFound = false;

    // Grab the list of seedling segments found in the graph connected to the junction coordinate of junctionKey
    std::pair<SeedlingSegmentGraph::iterator, SeedlingSegmentGraph::iterator> segments = segmentGraph->equal_range(junctionKey);
    SeedlingSegmentGraph::const_iterator segment;

    // Go through the list to see if the seedling segment in question exists in the graph
    for(segment = segments.first; segment != segments.second; ++segment) {
        // Do a simple pointer check to see if the two seedling segments are the same
        if(segment->second == seedling) {
            segmentFound = true;
            break;
        }
    }

    return segmentFound;
}

/* method getNeighbors
* Given x and y coordinates, populates the vector containing neighboring pixel coordinates.
* 
* @param x the x coordinate of the center pixel
* @param y the y coordinate of the center pixel
* @param neighbors vector of neighboring pixels that are part of the foreground 
*                  (e.g. connected)
*/
void UniformityAnalysis::getNeighbors(int x, int y, std::vector<coordinate> &neighbors) {
    // Temp pointer to current row in seedling image
    unsigned char *currentRow;

    // Check the 3x3 neighborhood around (x, y) in the seedlings image to discover unvisited neighbors
    for(int i = y + 1; i > y - 2; i--) {
        currentRow = (unsigned char *)(_seedlings->imageData + i * _seedlings->widthStep);
        for(int j = x - 1; j < x + 2; j++) {
            if(j == x && i == y) {
                continue;
            }
            if(currentRow[j] > 0 && isNotVisited(j, i)) {
                neighbors.push_back(std::make_pair(j, i));
            }
        }
    }
}

/* method segmentIncreasesLength
* Ensures that segment second increases the length of segment first by checking to see if
* it:  1) connects to the first segment, and 2) grows the segment upward
* 
* @param first   the seedling segment being grown
* @param second  the seedling segment being tested
* @return lengthIncreses  true if the segment meets the above criteria, false otherwise
*/
bool UniformityAnalysis::segmentIncreasesLength(SeedlingSegment *first, SeedlingSegment *second) {
    bool lengthIncreases = false;

    // Get the endpoints of both segments
    coordinate highestPointOfFirst = first->getMaxCoordinate();
    coordinate endOfFirst = first->getSegment().back();
    coordinate endOfSecond = second->getSegment().back();
    coordinate startOfSecond = second->getSegment().front();

    // Test to see of the second segment endpoint begins after the endpoint of the first.
    // Also check to make sure the segments are in fact connected.
    if( startOfSecond.first == endOfFirst.first && 
        startOfSecond.second == endOfFirst.second && 
        highestPointOfFirst.second > endOfSecond.second) {
            lengthIncreases = true;
    }

    return lengthIncreases;
}

/* method selectBestSegment
* Given the main seedling segment being grown and a vector of seedling segment choices,
* selects which of the choices to add to the main segment by picking the one which will
* have the smallest angle between the main segment and the candidate segments.
* 
* @param main     the main seedling segment
* @param choices  a vector of seedling segments connected to the main segment
*                 that will be evaluated using the above criteria
* @return bestCandidate  the seedling segment amongst the vector of choices which has
*                        the smallest angle between it and the main segment
*/
SeedlingSegment* UniformityAnalysis::selectBestSegment(CompleteSeedling *seedling, 
                                                       std::vector<SeedlingSegment*> &choices) {
    // Initialize local scratch and return variables
    double minAngle = 2 * M_PI, currentAngle;
    SeedlingSegment *bestCandidate = NULL;

    // Get the main segment from the seedling
    SeedlingSegment *main = seedling->getRootSegmentAt(-1);

    // If the main segment is growing north, but the segment choices are not,
    // use the subsegment of the main segment involving the second to last
    // segment
    bool useSubSegment = false;
    if(main->getEndDirection() == NORTH) {
       useSubSegment = true;
       for(unsigned int i = 0; i < choices.size() && useSubSegment; i++) {
           // Possibly change this to check both start and end directions are NORTH
           if(choices[i]->getStartDirection() == NORTH && choices[i]->getStartDirection() == NORTH) {
               useSubSegment = false;
           }
       }
    }

    // Select the appropriate segment for comparisons based on the previous check
    if(useSubSegment) {
       main = seedling->getRootSegmentAt(-2);
    }

    // Go through all seedling candidates and see which one can be connected to 
    // the main seedling segment with the minimum angle between both segments.
    for(unsigned int i = 0; i < choices.size(); i++) {
       currentAngle = measureAngle(main, choices[i]);
       if(currentAngle < minAngle) {
           minAngle = currentAngle;
           bestCandidate = choices[i];
       }
    }

    return bestCandidate;
}

/* method measureAngle
* Given two seeding segments, measures the internal crossing angle of their lines.
* 
* @param main  pointer to the main seedling segment
* @param candidate  pointer to the candidate seedling segment
* @return angle  double precision angle between the two segment lines in radians
*/
double UniformityAnalysis::measureAngle(SeedlingSegment *main, SeedlingSegment *candidate) {
    double angle = 2 * M_PI;

    // Compute indices for the start of the main segment and end of candidate segment
    int mainSize = main->getSegment().size(), candidateSize = candidate->getSegment().size();
    int main1index = mainSize > 20 ? mainSize - 10 : 0;
    int candidate2index = candidateSize > 20 ? 10 : candidateSize - 1;

    // Get the starting and ending coordinates of each seedling segment
    coordinate  main1 = main->getSegment().at(main1index),
        main2 = main->getSegment().back(),
        candidate1 = candidate->getSegment().front(),
        candidate2 = candidate->getSegment().at(candidate2index);

    // Compute the rise and run of each segment
    double  mainRise = main2.second - main1.second,
        mainRun  = main2.first - main1.first,
        candidateRise = candidate2.second - candidate1.second,
        candidateRun  = candidate2.first - candidate1.first;

    // Check both segments for infinite slope
    if(mainRun == 0 || candidateRun == 0) {
        // Both segments have an infinite slope
        if(mainRun == 0 && candidateRun == 0) {
            angle = 0.0;
        }
        // Main segment has infinite slope
        else if(mainRun == 0) {
            angle = atan(candidateRun / candidateRise);
        }
        // Candidate segment has infinite slope
        else {
            angle = atan(mainRun / mainRise);
        }
    }
    // Check both segments for zero slope
    else if(mainRise == 0 || candidateRise == 0) {
        // Both segments have a zero slope
        if(mainRise == 0 && candidateRise == 0) {
            angle = M_PI;
        }
        // Main segment has zero slope
        else if(mainRise == 0) {
            angle = atan(candidateRise / candidateRun);
        }
        // Candidate segment has zero slope
        else {
            angle = atan(mainRise / mainRun);
        }
    }
    // Both segments have non-zero and non-infinite slopes
    else {
        // Compute the slopes, their products and differences
        double mainSlope = mainRise / mainRun,
            candidateSlope = candidateRise/ candidateRun,
            product = mainSlope * candidateSlope,
            difference = mainSlope - candidateSlope;

        // Check for the special case where the product is -1
        if(product == -1) {
            angle = M_PI / 2;
        }
        /* Compute the slope using the formula:
        *
        *                     (mainSlope - candidateSlope)
        * angle = arctan( ------------------------------------ )
        *                   (1 - mainSlope * candidateSlope)
        */
        else {
            angle = atan( difference / (1 + product));
        }
    }

    // Convert the angle to degrees
    //angle *= (180 / M_PI);

    return abs(angle);
}

void UniformityAnalysis::pruneShortTerminalEdges(SeedlingSegmentGraph *segmentGraph) {
    SeedlingSegmentGraph::iterator segment;
    std::vector<SeedlingSegmentGraph::iterator> segmentsToDelete;
    int ep1Connections = 0, ep2Connections = 0;
    bool isInternalEdge = false;

    // Go through all graph edges (segments)
    for(segment = segmentGraph->begin(); segment != segmentGraph->end(); ++segment) {
        // Get the segment endpoints
        coordinate *ep1, *ep2;
        ep1 = &(segment->second->getSegment().back());
        ep2 = &(segment->second->getSegment().front());

        // Count the number of other edges endpoint 1 of this edge connects to
        _keyCreator.str(std::string());
        _keyCreator << ep1->first << "-" << ep1->second;
        ep1Connections = segmentGraph->count(_keyCreator.str()) - 1;

        // Count the number of other edges endpoint 2 of this edge connects to
        _keyCreator.str(std::string());
        _keyCreator << ep2->first << "-" << ep2->second;
        ep2Connections = segmentGraph->count(_keyCreator.str()) - 1;

        // Check the number of connections
        if(ep1Connections > 0 && ep2Connections > 0) {
            isInternalEdge = true;
        }
        else {
            isInternalEdge = false;
        }

        // If they meet the criteria for removal, add them to the removal list
        if( segment->second->getSegment().size() < MIN_ROOT_TIP_SIZE
            && segment->second->getMaxCoordinate().second > _seedlingBlobRect.y + 0.2 * _seedlingBlobRect.height
            && !isInternalEdge) {
                segmentsToDelete.push_back(segment);
        }
    }

    // Go through the removal list and remove marked edges from the graph
    for(size_t i = 0; i < segmentsToDelete.size(); i++) {
        segmentGraph->erase(segmentsToDelete[i]);
    }

    // Clear the key creator string stream as it may be used elsewhere
    _keyCreator.str(std::string());
}


/* method initVisitedPixelsArray
* Given a seedling blob, initializes an array the size of this blob's
* bounding rectangle and sets every element to zero.
* 
*/
void UniformityAnalysis::initVisitedPixelsArray(CBlob* blob) {
    // Get the bounding rectangle for this blob
    _seedlingBlobRect = blob->GetBoundingBox();

    // Make corrections for zero width or height blobs
    if(_seedlingBlobRect.width == 0) {
        _seedlingBlobRect.width = 1;
    }
    if(_seedlingBlobRect.height == 0) {
        _seedlingBlobRect.height = 1;
    }

    // Initialize the visited pixels array
    _visitedPixels = new bool[_seedlingBlobRect.width * _seedlingBlobRect.height];

    // Ensure it is completely zeroed out
    for(int y = 0; y < _seedlingBlobRect.height; y++) {
        for(int x = 0; x < _seedlingBlobRect.width; x++) {
            _visitedPixels[y * _seedlingBlobRect.width + x] = false;
        }
    }
}


/* method isNotVisited
* Given indices of pixels in the _seedlings image, isNotVisited checks the
* _visitedPixels array to determine if that pixel was previously processed
* 
* @param x  x index of the pixel in the_seedlings image
* @param y  y index of the pixel in the_seedlings image
* @return  true if the pixel has not been visited yet, false otherwise
*/
bool UniformityAnalysis::isNotVisited(int x, int y) {
    // Compute indices into visited pixels array
    int xIndex = x - _seedlingBlobRect.x;
    int yIndex = y - _seedlingBlobRect.y;
    int vpIndex = yIndex * _seedlingBlobRect.width + xIndex;

    // Return true if this pixel has not been visited, false otherwise
    if(0 <= vpIndex && vpIndex < _seedlingBlobRect.width * _seedlingBlobRect.height)
        return !_visitedPixels[vpIndex];
    else 
        return false;
}


/* method setVisited
* Given indices of pixels in the _seedlings image, setVisited
* sets the repsective _visitedPixels array value to true
* 
* @param x  x index of the pixel in the_seedlings image
* @param y  y index of the pixel in the_seedlings image
*/
void UniformityAnalysis::setVisited(int x, int y) {
    // Compute indices into visited pixels array
    int xIndex = x - _seedlingBlobRect.x;
    int yIndex = y - _seedlingBlobRect.y;
    int vpIndex = yIndex * _seedlingBlobRect.width + xIndex;

    // Set this pixel location as visited if it is within the bounds of the blob rectangle
    if(0 <= vpIndex && vpIndex < _seedlingBlobRect.width * _seedlingBlobRect.height)
        _visitedPixels[vpIndex] = true;
}


/* method findSeedlingStartPoint
* Using the bounds of the input seedling blob and the _seedlings image,
* find the start point of the next potential seedling within the image
* 
* @param seedling  the seedling blob
* @return  the coordinate from which to start building the seedling
*/
coordinate UniformityAnalysis::findSeedlingStartPoint(CBlob *seedling) {
    // Scan through the image looking for foreground
    // pixels since these pixels belong to seedlings.
    bool found = false;
    unsigned char *currentRow;
    coordinate startCoordinate(-1, -1);
    int widthStep = _seedlings->widthStep;
    unsigned char *imageData = (unsigned char *) _seedlings->imageData;

    // Search for the start of the seedling within the seedling image
    for(int y = (int)seedling->MaxY(); y >= (int)seedling->MinY() && !found; y--) {
        currentRow = (unsigned char *)(imageData + y * widthStep);

        // Search within the row
        for(int x = (int)seedling->MinX(); x <= (int)seedling->MaxX(); x++) {

            // Starting point has been found
            if(currentRow[x] > 0 && isNotVisited(x, y)) {
                startCoordinate.first = x;
                startCoordinate.second = y;
                found = true;
                break;
            }
        }
    }

    // Return the coordinate
    return startCoordinate;
}

/*
int UniformityAnalysis::roundToNearest5MM(int length) {
// Extract the ones digit from the length
int onesDigit = length % 10;

// Subtract the ones digit from the length to prepare for rounding
int tempLength = length - onesDigit;

// Test to see whether to round ones digit down to 0 or up to 5 mm
if(0 <= onesDigit && onesDigit < 5) {
// Check ones digit distance from 0 and 5
if(onesDigit > (5 - onesDigit)) {			
// (onesDigit - 0) distance > (5 - onesDigit) distance, round up to 5 mm
tempLength += 5;
}
}
// One's digit is greather than 5, test to see if tempLength should
// be rounded up to the next nearest 10 or left at the nearest 5 mm
else {
// Start at the nearest 5 mm before test to see if tempLength 
// should stay at nearest 5 mm or rounded up to nearest 10 mm
tempLength += 5;

// Check ones digit distance from 5 and 10
if((onesDigit - 5) > (10 - onesDigit)) {			
// (onesDigit - 5) distance > (10 - onesDigit) distance, round up to nearest 10 mm
tempLength += 5;
}
}

// Return the final rounded length
return tempLength;
}
*/


/* method computeUniformity
* Calculates the uniformity rating using the the following approach:
*   - Finds the _numGroups largest contiguous groups within _groups.
*   - Sums up their counts and returns this sum divided by the total
*     number of seedlings detected within the image.
* @return uniformity index computed as described above
*/
double UniformityAnalysis::computeUniformity(void) {
    double uniformity = 0.0;

    // Create a list containing (length, count) pairs sorted by group length
    std::vector<int> sortedGroups(_groupKeys.begin(), _groupKeys.end());
    std::sort(sortedGroups.begin(), sortedGroups.end());

    // Check to see if any seedlings should be combined together, this
    // is done to fix seedlings that are broken due to thresholding


    // Compare counts of each _numGroups contiguous groups to compute 
    // the uniformity. The final uniformity is the group count sums
    // divided by the total number of seedlings discovered;
    std::cout << "Number of Seedlings Found: " << _seedlingsFound << endl;
    unsigned int currentSum = 0, maxSum = 0;
    for(unsigned int i = 0; i < sortedGroups.size() - _numGroups; i++) {
        currentSum = 0;
        for(unsigned int j = i; j < i + _numGroups; j++) {
            currentSum += _groups.count(sortedGroups.at(j));
        }
        if(currentSum > maxSum) {
            maxSum = currentSum;
        }
    }

    // Compute the uniformity percentage as:
    // ( max group count sum / total seedling count ) * 100
    uniformity = ( (double) maxSum / (double) _seedlingsFound ) * 100;

    return uniformity;
}



/* method skeletonize
* Implements the two-pass Zhang and Suen skeletonization algorithm  as
* seen in Section 9.6 page 300 of Advanced Biomedical Image Analysis.
* The first pass marks pixels for deletion, while the second pass goes
* through and deletes the pixels.
*
* @param img The single channel image that should be skeletonized
* @return skeletonized version of input image
*/
IplImage* UniformityAnalysis::skeletonizeZS(IplImage *img) {
    // Initialize the iteration and pixels removed counters
    //CvScalar pixelToCheck, deletionValue = cvScalarAll(0.);
    int iteration = 0, numPixelsRemoved, widthStep = img->widthStep, x, y, bitCode;
    unsigned char *currentRow, *rowAbove, *rowBelow;

    // Grab the image width and height for later use
    int width = img->width, height = img->height;

    // Initialize the skeleton image
    IplImage *skeleton = cvCreateImage(cvSize(width, height), img->depth, img->nChannels);

    // Ensure a single channel image was provided
    if(img->nChannels != 1) {
        return skeleton;
    }

    // Initialize the array F, which is used to mark pixels for deletion
    int *F = new int[widthStep*height];

    // Copy the input image data into the skeleton image
    cvCopy(img, skeleton, 0);

    // Begin iteratively marking and deleting pixels until no changes occur
    do {
        // Reset the pixels removed counter
        numPixelsRemoved = 0;
        // Iterate over the entire image, marking pixels for deletion
        // Iterate over the rows
        for(y = 1; y < height - 1; y++) {
            // Pointer to the start of the current row, row above and row below
            currentRow = (unsigned char *)(skeleton->imageData + y * skeleton->widthStep);
            rowAbove = (unsigned char *)(skeleton->imageData + (y - 1) * skeleton->widthStep);
            rowBelow = (unsigned char *)(skeleton->imageData + (y + 1) * skeleton->widthStep);

            // Iterate over the columns in the current row
            for(x = 1; x < width - 1; x++) {
                // Initially, the pixel should not be marked for deletion
                F[y * widthStep + x] = 0; 

                // Only perform further processing if the current pixel is in the foreground
                if(currentRow[x] > 0) {
                    // Reset the bit code used to determine 
                    // if the center pixel should be removed
                    bitCode = 0;

                    // Use the 8-connected neighborhood to compute the bit code
                    if(rowAbove[x] > 0)			bitCode += 1;
                    if(rowAbove[x + 1] > 0)		bitCode += 2;
                    if(currentRow[x + 1] > 0)	bitCode += 4;
                    if(rowBelow[x + 1] > 0)		bitCode += 8;
                    if(rowBelow[x] > 0)			bitCode += 16;
                    if(rowBelow[x - 1] > 0)		bitCode += 32;
                    if(currentRow[x - 1] > 0)	bitCode += 64;
                    if(rowAbove[x - 1] > 0)		bitCode += 128;

                    // Determine if the center pixel should be marked for deletion
                    if( (fateTable[bitCode] == 3) ||
                        (fateTable[bitCode] == 2 && ((iteration & 1) == 0)) ||
                        (fateTable[bitCode] == 1 && ((iteration & 1) == 1)) )
                        F[y * widthStep + x] = 1;
                }
            }
        }

        // Iterate over the image now and delete the marked pixels
        // Iterate over the rows
        for(y = 1; y < height - 1; y++) {
            // Pointer to the start of the current row
            currentRow = (unsigned char *)(skeleton->imageData + y * img->widthStep);
            // Iterate over the columns in the current row
            for(x = 1; x < width - 1; x++) {
                // If the current pixel was marked for deletion, delete it
                if(F[y * widthStep + x] == 1) {
                    // Increment the pixels removed counter
                    numPixelsRemoved++;

                    // Delete the pixel
                    currentRow[x] = 0;
                }
            }
        }

        // Increment the iteration counter
        iteration++;
    } while(numPixelsRemoved > 0);

    // Free F matrix memory
    delete [] F;

    // Return the skeletonized image
    return skeleton;
}



void UniformityAnalysis::thresholdRange(IplImage *src, IplImage *dst, int lowerThresh, int upperThresh, int maxValue, bool keepSrcValue) {
    assert(src->width == dst->width && dst->height == src->height);
    unsigned char *currentRowSrc, *currentRowDst;
    for(int y = 0; y < src->height; y++) {
        // Pointer to the start of the current row
        currentRowSrc = (unsigned char *)(src->imageData + y * src->widthStep);
        currentRowDst = (unsigned char *)(dst->imageData + y * dst->widthStep);

        // Iterate over the columns in the current row
        for(int x = 0; x < src->width; x++) {
            if(lowerThresh <= currentRowSrc[x] && currentRowSrc[x] <= upperThresh) {
                if(keepSrcValue) {
                    currentRowDst[x] = currentRowSrc[x];
                }
                else {
                    currentRowDst[x] = maxValue;
                }
            }
            else {
                currentRowDst[x] = 0;
            }
        }
    }
}



void UniformityAnalysis::invertImage(IplImage *src, IplImage *dst) {	
    assert(src->width == dst->width && dst->height == src->height);
    assert(src->nChannels == 1 && dst->nChannels == 1);
    int x, y;
    unsigned char *currentRowSrc, *currentRowDst;
    for(y = 0; y < src->height; y++) {
        // Pointer to the start of the current row
        currentRowSrc = (unsigned char *)(src->imageData + y * src->widthStep);
        currentRowDst = (unsigned char *)(dst->imageData + y * dst->widthStep);

        // Iterate over the columns in the current row
        for(x = 0; x < src->width; x++) {
            currentRowDst[x] = 255 - currentRowSrc[x];
        }
    }
}


void UniformityAnalysis::bitwiseAnd(IplImage *src1, IplImage *src2, IplImage *dst) {
    assert(src1->width == src2->width && src2->height == src1->height);
    assert(src1->width == dst->width && dst->height == src1->height);
    assert(src1->nChannels == 1 && dst->nChannels == 1 && src2->nChannels == 1);
    int x, y;
    unsigned char *currentRowSrc1, *currentRowSrc2, *currentRowDst;
    for(y = 0; y < src1->height; y++) {
        // Pointer to the start of the current row
        currentRowSrc1 = (unsigned char *)(src1->imageData + y * src1->widthStep);
        currentRowSrc2 = (unsigned char *)(src2->imageData + y * src2->widthStep);
        currentRowDst = (unsigned char *)(dst->imageData + y * dst->widthStep);

        // Iterate over the columns in the current row
        for(x = 0; x < src1->width; x++) {
            currentRowDst[x] = currentRowSrc1[x] & currentRowSrc2[x];
        }
    }
}

void UniformityAnalysis::printSeedlingInfo(CompleteSeedlings *seedlings) {
    for(unsigned int i = 0; i < seedlings->size(); i++) {
        std::cout << "-----------------------------------" << std::endl;
        std::cout << "   Printing Info For Seedling #" << i << std::endl;
        std::cout << "-----------------------------------" << std::endl;
        seedlings->at(i)->printSeedlingInfo();
    }
}

/******************************************************
* Experimental second approach to building seedlings *
******************************************************/

SeedlingSegmentVector *UniformityAnalysis::findRootTips(SegmentVectorGraphPair *svg) {
    // Initialize root tip seedling segment container
    SeedlingSegmentVector *rootTipCandidates = new SeedlingSegmentVector();

    // Grab the segment vector out of the svg pair
    SeedlingSegmentVector *segmentVector = svg->first;

    // Go through the segment vector and determine the root tip candidates
    for(size_t i = 0; i < segmentVector->size(); i++) {
        // Grab the segment's coordinates list 
        std::vector<coordinate> *coords = &(segmentVector->at(i)->getSegment());

        // Check that it meets the minimum requirements for a root
        // tip before checking other more complicated requirements
        if(coords->size() > MIN_ROOT_TIP_SIZE) {
            // Determine which endpoint, if any, of this segment is the terminal one
            coordinate *terminalEndpoint = NULL;

            // Segment ending pixel is the terminal endpoint
            if(isJunctionPixel(coords->front())) {
                terminalEndpoint = &(coords->back());
            }
            // Segment starting pixel is the terminal endpoint
            else if(isJunctionPixel(coords->back())) {
                terminalEndpoint = &(coords->front());
            }

            // There is a chance that neither endpoint is a terminal pixel,
            // which means that this segment is a connecting segment and so
            // cannot be a root tip. Exclude it from any further processing.
            if(terminalEndpoint != NULL) {
                // The root tip is generally found in the bottom 2/3rds of the seedling.
                // This stage is to ensure that whenever there are overlapping segments,
                // all potential seedlings are found, while ignoring segments that can't
                // be root tips.
                if((_seedlingBlobRect.y + 0.667 * _seedlingBlobRect.height) <= terminalEndpoint->second) {
                    // This segment passed all requirements and is, therefore, a root tip candidate
                    rootTipCandidates->push_back(segmentVector->at(i));
                }
            }
        }
    }

    return rootTipCandidates;
}

CompleteSeedlings* UniformityAnalysis::createSeedlings(SegmentVectorGraphPair *svg) {
    // Initialize the group container for all seedlings
    CompleteSeedlings *seedlingGroups = new CompleteSeedlings();	

    // Prior to sorting the seedling segment vector by Y coordinates, make sure that all seedlings are growing
    // in the same direction. Thus, if a seedling's starting Y coordinate is less than its ending Y coordinate,
    // then reverse the seedling segment and make the old end coordinate the new start coordinate.
    correctSeedlingEndpoints(svg->first);

    // Check the special case where there is only one seedling segment
    if(svg->first->size() == 1) {
        seedlingGroups->push_back(new CompleteSeedling());
        seedlingGroups->front()->addToPrimaryRoot(svg->first->front());
    }
    // Otherwise run through the normal processing steps
    else {
        // Get the segment graph and root tip candidate vector from the svg pair
        SeedlingSegmentVector *rootTipCandidates = findRootTips(svg);
        SeedlingSegmentGraph *segmentGraph = svg->second;

        // Begin creating the groups
        for(size_t i = 0; i < rootTipCandidates->size(); i++) {
            // Create a group for this root candidate
            seedlingGroups->push_back(new CompleteSeedling());
            seedlingGroups->at(i)->addToPrimaryRoot(rootTipCandidates->at(i));
            createGroup(segmentGraph, seedlingGroups->at(i));
        }
    }

    // return the seedling groups for further processing
    return seedlingGroups;
}

void UniformityAnalysis::createGroup(SeedlingSegmentGraph *segmentGraph, CompleteSeedling *seedling) {
    // Pointer to the next segment to examine
    SeedlingSegment* mainSegment = seedling->getRootSegmentAt(0);

    // Keep a list of the junctions that should be/have been explored while building the seedling
    std::list<coordinate> junctions;
    std::vector<coordinate> visitedJunctions;

    // Add the first junction to be examined
    junctions.push_back(mainSegment->getSegment().back());

    // Build a vector of all segments connected to this junction, if they exist
    while(!junctions.empty()) {
        // Get the next junction to check
        coordinate nextJunction = junctions.back();
        junctions.pop_back();
        visitedJunctions.push_back(nextJunction);

        // Build the key string
        _keyCreator << nextJunction.first << "-" << nextJunction.second;

        // Check all connected segments and grow the seedling with the best one
        addSegmentsToGroup(segmentGraph, seedling, _keyCreator.str(), junctions);

        // Clear the key creator string stream
        _keyCreator.str("");

        // Update the main segment and add the next junction pixel that should be checked, if any
        mainSegment = seedling->getRootSegmentAt(-1);
        coordinate endOfMainSegment = mainSegment->getSegment().back();
        bool addJunction = true;

        // Search through the visited junctions to ensure we don't look at the same junction twice
        for(unsigned int i = 0; i < visitedJunctions.size(); i++) {
            if( visitedJunctions[i].first == endOfMainSegment.first &&
                visitedJunctions[i].second == endOfMainSegment.second) {
                    addJunction = false;
                    break;
            }
        }

        // Add the junction if we havent checked it already
        if(addJunction) {
            junctions.push_back(mainSegment->getSegment().back());
        }

        // Update the main segment
        mainSegment = seedling->getRootSegmentAt(-1);
    }
}

void UniformityAnalysis::addSegmentsToGroup(SeedlingSegmentGraph* segmentGraph, CompleteSeedling* seedling, std::string &mainSegmentKey, std::list<coordinate> &junctions) {

    // Graph iterator used to go through all seedling segment edges
    SeedlingSegmentGraph::const_iterator currentSegment;

    // Get the main segment from the connectedSegments vector
    SeedlingSegment* mainSegment = seedling->getRootSegmentAt(-1);

    // A pair of iterators used to obtain pointers to the seedling
    // segments connected to the current junction key
    SeedlingSegmentEdges edges = segmentGraph->equal_range(mainSegmentKey);

    // Go through all the segments connected to the current junction
    for(currentSegment = edges.first; currentSegment != edges.second; ++currentSegment) {
        // Clear the keyCreator for use with the current segment
        _keyCreator.str(std::string());

        // Skip the segment connectiong to this junction that is already part of the seedling
        if(seedling->isPartOfSeedling(currentSegment->second)) {
            continue;
        }

        // Get the endpoint of the segment that is not connect to this junction
        coordinate endpoint = currentSegment->second->getSegment().back();

        // Create the segment graph junction key for this endpoint
        _keyCreator << endpoint.first << "-" << endpoint.second;

        // Ensure the key is not the same as the key mainSegmentKey
        if(mainSegmentKey == _keyCreator.str()) {
            _keyCreator.str(std::string());
            endpoint = currentSegment->second->getSegment().front();
            _keyCreator << endpoint.first << "-" << endpoint.second;
        }

        // Consider the segment if it increases the length of the main segment
        if(segmentIncreasesLength(mainSegment, currentSegment->second)) {
            // If it is a terminal endpoint and in the top third
            // of the seedling blob add it to the seedling group
            if(segmentGraph->count(_keyCreator.str()) < 2 &&
               currentSegment->second->getMaxCoordinate().second < (_seedlingBlobRect.y + 0.33 * _seedlingBlobRect.height)) {
                    seedling->addToPrimaryRoot(currentSegment->second);
            }
            else if(currentSegment->second->getSegment().size() > MIN_ROOT_TIP_SIZE) {
                seedling->addToPrimaryRoot(currentSegment->second);
            }

            // Add the endpoint to the junctions queue
            junctions.push_back(endpoint);
        }
        // If it is a connecting segment without inflection points
        else if(segmentGraph->count(_keyCreator.str()) > 1 && !currentSegment->second->hasInflectionPoint()) {
            seedling->addToPrimaryRoot(currentSegment->second);

            // Add the endpoint to the junctions queue
            junctions.push_back(endpoint);
        }
        // This segment is most likely a secondary root segment. However, before
        // adding it as a secondary segment, check that it is a terminal segment.
        else if(currentSegment->second->getSegment().size() < MIN_ROOT_TIP_SIZE){
            seedling->addAsSecondaryRoot(currentSegment->second);
        }
    }
}


const int UniformityAnalysis::fateTable[256] = {
    0,0,0,1,0,3,1,1,0,0,0,0,2,2,3,3,0,0,0,0,3,0,0,0,2,0,0,0,2,0,3,2,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,2,0,2,0,3,0,0,0,3,0,3,2,
    0,3,0,1,0,0,0,0,0,0,0,0,0,0,0,0,3,0,0,0,0,0,0,0,2,0,0,0,0,0,0,0,
    2,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,0,0,0,0,3,0,0,0,2,0,2,0,
    0,1,0,3,0,1,0,3,0,0,0,0,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    1,1,0,3,0,0,0,1,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    3,3,0,3,0,0,0,1,0,0,0,0,0,0,0,0,3,1,0,1,0,0,0,0,3,1,0,0,2,0,0,0
};

const CvScalar UniformityAnalysis::colors[20]  = {
    cvScalar(0xff, 0xb3, 0x00), cvScalar(0x80, 0x3e, 0x75),
    cvScalar(0xff, 0x68, 0x00), cvScalar(0xa6, 0xbd, 0xd7),
    cvScalar(0xc1, 0x00, 0x20), cvScalar(0xce, 0xa2, 0x62),
    cvScalar(0x81, 0x70, 0x66), cvScalar(0x00, 0x7d, 0x34),
    cvScalar(0xf6, 0x76, 0x8e), cvScalar(0x00, 0x53, 0x8a),
    cvScalar(0xff, 0x7a, 0x5c), cvScalar(0x53, 0x37, 0x7A),
    cvScalar(0xFF, 0x8E, 0x00), cvScalar(0xB3, 0x28, 0x51),
    cvScalar(0xF4, 0xC8, 0x00), cvScalar(0x7F, 0x18, 0x0D),
    cvScalar(0x93, 0xAA, 0x00), cvScalar(0x59, 0x33, 0x15),
    cvScalar(0xF1, 0x3A, 0x13), cvScalar(0x23, 0x2C, 0x16)
};

const int UniformityAnalysis::validJunctionSums[24] = {
    21, 37, 41, 42,57, 69, 73, 74, 81, 82, 83, 84, 85, 
    138, 146, 147, 162, 164, 165, 168, 169, 170, 172, 228
};

const set<int> UniformityAnalysis::junctionSumSet(validJunctionSums, validJunctionSums + 22);